par(mfrow=c(2, 2))
data1_14 <- read.csv("data/时间序列分析——基于R（第2版）案例数据/csv/A1_14.csv")
ts_GNP <- ts(data1_14$GNP)

# 时序图
plot(ts_GNP, type="o", pch=5, col='#39CBB4')


# ADF检验
# install.packages("aTSA")
library(aTSA)
adf.test(ts_GNP, nlag=3)

# 一阶差分
dif_GNP<-diff(ts_GNP)
plot(dif_GNP, type="o", pch=5, col='#39CBB4')

# ADF检验
adf.test(dif_GNP, nlag=3)


# 白噪声，纯随机检验
for( k in 1:3) print(Box.test(dif_GNP, lag=6*k, type="Ljung-Box"))


# 自相关图
acf(dif_GNP, lag.max=30)
pacf(dif_GNP, lag.max=30)

# 拟合 ARMA(1, 1) 模型
model11 <- arima(dif_GNP, order=c(1, 0, 0), method="ML")
model12 <- arima(ts_GNP, order=c(1, 1, 0), method="ML")

# 或者 ARMA(1, 1) 模型
# install.packages("forecast")
library(forecast)
model21 <- Arima(dif_GNP, order=c(1, 0, 0), method="ML", include.drift = T)
model22 <- Arima(ts_GNP, order=c(1, 1, 0), method="ML", include.drift = TRUE)


model12
model22

# 预测未来10期
forecast(model12, h = 10)
forecast(model22, h = 10)

